home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / hamradio / sgp4_pl2.zip / SOLAR.PAS < prev    next >
Pascal/Delphi Source File  |  1992-10-01  |  3KB  |  94 lines
  1. Unit Solar;
  2. {           Author:  Dr TS Kelso }
  3. { Original Version:  1990 Jul 29 }
  4. { Current Revision:  1992 Sep 30 }
  5. {          Version:  1.20 }
  6. {        Copyright:  1990-1992, All Rights Reserved }
  7. {$N+}
  8.  
  9. INTERFACE
  10.   Uses SGP_Math;
  11.  
  12. const
  13.   eclipsed : boolean = false;
  14.   show_vis : boolean = false;
  15.  
  16. var
  17.   civil,
  18.   nautical,
  19.   astronomical : double;  {Twilight elevations}
  20.   solar_pos    : vector;
  21.  
  22. Procedure Calculate_Solar_Position(time : double;
  23.                        var solar_vector : vector);
  24. Function Depth_of_Eclipse(time : double;
  25.                             r1 : vector) : double;
  26.  
  27. IMPLEMENTATION
  28.   Uses SGP_Intf,SGP_Time;
  29.  
  30. const
  31.   sr = 696000.0;       {Solar radius - kilometers (IAU 76)}
  32.   AU = 1.49597870E8;   {Astronomical unit - kilometers (IAU 76)}
  33.  
  34. Procedure Calculate_Solar_Position(time : double;
  35.                        var solar_vector : vector);
  36.   var
  37.     mjd,year,T,M,L,e,C,O,Lsa,nu,R,eps : double;
  38.     ob                                : vector;
  39.   begin
  40.   mjd := time - 2415020.0;
  41.   year := 1900 + mjd/365.25;
  42.   T := (mjd + Delta_ET(year)/secday)/36525.0;
  43.   M := Radians(Modulus(358.47583 + Modulus(35999.04975*T,360.0)
  44.      - (0.000150 + 0.0000033*T)*Sqr(T),360.0));
  45.   L := Radians(Modulus(279.69668 + Modulus(36000.76892*T,360.0)
  46.      + 0.0003025*Sqr(T),360.0));
  47.   e := 0.01675104 - (0.0000418 + 0.000000126*T)*T;
  48.   C := Radians((1.919460 - (0.004789 + 0.000014*T)*T)*Sin(M)
  49.      + (0.020094 - 0.000100*T)*Sin(2*M) + 0.000293*Sin(3*M));
  50.   O := Radians(Modulus(259.18 - 1934.142*T,360.0));
  51.   Lsa := Modulus(L + C - Radians(0.00569 - 0.00479*Sin(O)),twopi);
  52.   nu := Modulus(M + C,twopi);
  53.   R := 1.0000002*(1 - Sqr(e))/(1 + e*Cos(nu));
  54.   eps := Radians(23.452294 - (0.0130125 + (0.00000164 - 0.000000503*T)*T)*T
  55.        + 0.00256*Cos(O));
  56.   R := AU*R;
  57.   solar_vector[1] := R*Cos(Lsa);
  58.   solar_vector[2] := R*Sin(Lsa)*Cos(eps);
  59.   solar_vector[3] := R*Sin(Lsa)*Sin(eps);
  60.   solar_vector[4] := R;
  61.   end; {Procedure Calculate_Solar_Position}
  62.  
  63. Function Depth_of_Eclipse(time : double;
  64.                             r1 : vector) : double;
  65.   var
  66.     r1_r1,r1_r2,r2_r2,k,d,ds : double;
  67.     r2                       : vector;
  68.   begin
  69.   Magnitude(r1);
  70.   Calculate_Solar_Position(time,r2);
  71.   solar_pos := r2;
  72.   r1_r1 := Sqr(r1[4]);
  73.   r1_r2 := -Dot(r1,r2);
  74.   r2_r2 := Sqr(r2[4]);
  75.   k := r1_r2/r2_r2;
  76. {Calculate perpendicular distance from anti-solar vector}
  77.   d := Sqrt(r1_r1 - Sqr(r1_r2)/r2_r2);
  78. {Calculate shadow distance ds}
  79.   ds := xkmper + k * (sr - xkmper);
  80. {If d < ds, then satellite is in eclipse}
  81.   if (k > 0) and (d < ds) then
  82.     eclipsed := true
  83.   else
  84.     eclipsed := false;
  85.   Depth_of_Eclipse := d - ds
  86.   end; {Function Depth_of_Eclipse}
  87.  
  88.   begin
  89.   civil        := Radians(-6);
  90.   nautical     := Radians(-12);
  91.   astronomical := Radians(-18);
  92.  
  93. end.
  94.